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ABSTRACT 

The aim of this work is to compare the orbital dynamics in three different mod¬ 
els describing the properties of a star cluster rotating around its parent galaxy in 
a circular orbit. In particular, we use the isochrone and the Hernquist potentials to 
model the spherically symmetric star cluster and we compare our results with the 
corresponding ones of a previous work in which the Plummer model was applied for 
the same purpose. Our analysis takes place both in the configuration (x,y) and in 
the phase (x, x) space in order to elucidate the escape process as well as the overall 
orbital properties of the tidally limited star cluster. We restrict our investigation into 
two dimensions and we conduct a thorough numerical analysis distinguishing between 
ordered and chaotic orbits as well as between trapped and escaping orbits, considering 
only unbounded motion for several energy levels above the critical escape energy. It 
is of particular interest to determine the escape basins towards the two exit channels 
(near the Lagrangian points L\ and L^) and relate them with the corresponding escape 
times of the orbits. 
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1 INTRODUCTION 

We consider the simplest case of a star cluster in a circular 
orbit around its parent galaxy. Moreover, the star cluster is 
embedded in the tidal filed of the parent galaxy with lin¬ 
ear tidal forces. Usually we use the term “tidally limited” 
star cluster which means, roughly speaking, that the tidal 
field imposes a boundary, outside of which the mass density 
vanishes, while the star cluster is being captured within the 
tidal limit. Stars escape from tidally limited star clusters 
following a two stages process. During the first stage, stars 
are being moved by certain physical processes into the es¬ 
caping phase space, while in the second stage they escape 
through the channels in the open equipotential surface. The 
time needed for a star to complete stage two is mainly re¬ 
lated with the energy of the star. Fukushige & Heggie (2000) 
derived an expression for the escape times of stars in second 
stage which have just finished stage one, while Baumgardt 
(2001) exploited these results in order to address the im¬ 
portant issue of dissolution time of star clusters moving in 
circular orbits around their parent galaxy. 

Escaping particles from dynamical systems is a 
subject to which has been devoted many studies over 
the years. Especially the issue of escapes in Hamil- 
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tonian systems is directly related to the problem of 
chaotic scattering which has been an active field of 
research over the last decades and it still remains 
open (e.g., Bleher et al. 1988, 1989; Jung h Scholz 

1988; Contopoulos h Kaufmann 1992; Motter h Lai 
2002; Seoane et al. 2006, 2007; Seoane & Sanjuan 2008; 
Seoane et al. 2009; Seoane h Sanjuan 2010). The problem 
of escape is a classical problem in simple Hamiltonian non¬ 
linear systems (e.g., Aguirre et al. 2001; Aguirre & Sanjuan 
2003; Aguirre et al. 2009; Barrio et al. 2008; Blesa et al. 
2012; Zotos 2014, 2015c) as well as in dynamical as¬ 
tronomy (e.g., Hut & Bahcall 1983; Benet et al. 1996, 
1998; de Moura h Letelier 2000; Nagler 2004, 2005; 

de Assis & Terra 2014; Zotos 2012, 2015b, d). 

The chaotic dynamics of a star cluster embedded in the 
tidal field of a parent galaxy was investigated in Ernst et al. 
(2008) (hereafter Paper I). Conducting a thorough scanning 
of the available phase space the authors managed to obtain 
the basins of escape and the respective escape rates of the 
orbits, revealing that the higher escape times correspond 
to initial conditions of orbits near the fractal basin bound¬ 
aries. The investigation was expanded into three dimensions 
in Zotos (2015a) where we revealed the escape mechanism 
of three-dimensional orbits in a tidally limited star cluster. 
Furthermore, Ernst & Peters (2014) explored the escape dy¬ 
namics in the close vicinity of and within the critical area 
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in a two-dimensional barred galaxy potential, identifying the 
escape basins both in the phase and the configuration space. 

In Paper I the authors used the Plummer potential in 
order to model the properties of the star cluster. However 
there are also other spherically symmetric potentials which 
can also be used for the same purpose. For example there is 
the Henon’s isochrone and the Hernquist potential, where 
both of them have been previously used to describe the 
properties of a star cluster (e.g., Brasser et al. 2012; Adams 
2006). In the present work, we shall use these two potentials 
in an attempt to compare the orbital dynamics of the star 
cluster with the corresponding outcomes of Paper I. 

The paper is organized as follows: A detailed presenta¬ 
tion of the tidal approximation model is included in Section 
2. In Section 3 we describe all the computational techniques 
we used in order to determine the nature (bounded or un¬ 
bounded) of the orbits. The following Section contains all 
the numerical results showing how the value of the total 
orbital energy influences the orbital content of the models. 
Our paper ends with Section 5, where the conclusions and 
the discussion of this work are given. 


2 THE TIDAL APPROXIMATION MODEL 


Let us briefly recall the properties of the tidal approxima¬ 
tion theory. This model allow us to assume that the star 
cluster moves with constant angular velocity uj around the 
center of the parent galaxy. Using this model we can apply 
epicyclic approximation to calculate linear tidal forces act¬ 
ing on all members of the cluster. Undoubtedly, the most 
appropriate system of coordinates is an accelerated rotating 
reference frame (Chandrasekhar 1942) in which both the the 
star cluster and the galactic center are at rest, while the ori¬ 
gin of the coordinates is at the center of the star cluster 
which is located at the extremum of the effective potential 
of the galaxy, so that the position of the galactic center is 
(—R g ,0,0), where R g is the radius of the circular orbit. If 
the motion of the star cluster around its parent galaxy is 
more general (i.e., following an elliptical orbit), then no in¬ 
tegral of motion comparable to the Jacobi integral is known 
and therefore all the numerical computations become more 
difficult to be performed. 

In order to model the properties of the star cluster we 
use two types of spherically symmetric potentials: 

(a). The Henon’s isochrone potential (Henon 1959) 




GM c i 

ri + s/r 1 +r t 2 ’ 


( 1 ) 


where r 2 = x 2 + y 2 + z 2 , G is the gravitational constant, M c \ 
is the total mass of the star cluster, while n is the Isochrone 
radius similar to the Plummer radius. 

(b). The Hernquist potential (Hernquist 1990) 


4> h (x,y,z) 


GM C \ 
ru + r’ 


( 2 ) 


where m is the Hernquist radius. Potential (2) is in fact a 
Dehnen model (Dehnen 1993) with 7=1. 

The total effective potential 

<h e ff (x, y , z) = $ c i(x, y, z) + i (ft 2 - 4a; 2 ) x 2 + ^z/ 2 z 2 , (3) 


where <f> c i is either the isochrone or the Hernquist potential. 


A fundamental scale length of the star cluster is the 
tidal radius (King 1962) which is defined as 


( GMjr t ) \ 

\4uj 2 -k 2 J ’ 


( 4 ) 


where M{r t ) is the mass contained within the tidal radius 
and usually is smaller than the total mass of the star cluster 
M c 1 . 

We use the dimensionless system of units introduced in 
Paper I according to which G = uj — rt — 1. The formula 
which relates the physical units of time, mass and radius is 
G = 1/222.3 pc 3 /M 0 /Myr 2 . The characteristic frequencies 
uj, ft and v on the other hand, are related to the galactic 
gravitational potential and in the solar neighborhood can 
be expressed as functions of the Oort’s constants (see e.g., 
Binney & Tremaine 2008) as follows 

u 2 = (A- Bf , (5) 

k 2 = -4 B ( A - B) 2 , ( 6 ) 

v 2 = 4tt Gp s + 2 (A 2 - B 2 ) , (7) 


where p g is the local value of the galactic density 
(Holmberg & Flynn 2000). Exploiting the numerical val¬ 
ues of the Oort’s constants provided in Feast & Whitelock 
(1997), we obtain the dimensionless quantities ft 2 /uj 2 ~ 1.8 
and v 2 /uj 2 ~ 7.6. 

In both models we choose that M c \ = 2.7, while the 
mass contained within the tidal radius is M(r*t) = 2.2, which 
means that about 18.5% of the total mass lies beyond the 
tidal radius. For the value M(r t ) = 2.2 we obtain form Eq. 
(4) that the tidal radius is r t = 1 since 4a; 2 — ft 2 = 2.2. 
The second step is to set the values of the isochrone and 
the Hernquist radius. For this task we need a formula that 
relates the mass with the scale radius. We know that in 
spherical potentials 

M GM(r) 

dr r 2 ’ ^ 


so the cumulative mass distribution is 


M(r) 


G dr ' 


Inserting Eqs. (1) and (2) in formula (9) we have 


( 9 ) 


M(r) 


_ M c ir 3 _ 

y^TT 2 (n + yr 2 + r 2 ) 


( 10 ) 


for the isochrone potential and 


M(r ) = 


M c ir 2 

(r+ r H ) 2 ' 


( 11 ) 


for the Hernquist potential (see also Eq. (3) in Hernquist 
(1990)). We already know that M c \ = 2.7, while for r = 
rt = 1 it is M(r = r t ) — 2.2. Replacing these values in 
Eqs. (10) and (11) and solving for r while keeping always 
the positive solution we finally obtain that n — 0.10, while 
m = 0.107823. These values secure that the two Lagrangian 
saddle points L\ and L 2 are located exactly at ( x,y,z ) = 
( — 1,0,0) and ( x,y,z ) = (1,0,0), respectively. 

Taking into account the relationships connecting the 
tidal forces with the epicyclic frequency k and the vertical 
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Figure 1. Evolution of the mass density of the three spherically 
symmetric potentials (Plummer, isochrone and Hernquist) as a 
function of the radius r. 


frequency v the basic equations of motion are 


x — — 


y = - 


z — — 


d£ci 

dx 

d£ci 

dy 

<9$d 

dz 


— (k 2 — 4w 2 ) x + 2 ujy, 

— 2ljx, 


( 12 ) 

(13) 

(14) 


where, as usual, the dot indicates derivative with respect 
to the time. We observe that both centrifugal and Coriolis 
forces appear in the above equations of motion due to the 
fact that the coordinate system is rotating. 

In addition, the equations of motion admit the following 
isolating integral of motion 


C = I (x 2 + y 2 + z 2 ) + 3> e ff (*, y, z) = E, 


(15) 


where x, y and z are the momenta per unit mass, conjugate 
to x, y and z, respectively, while E is the numerical value 
of the Jacobian integral, which is conserved. The numerical 
value of the effective potential at the two Lagrangian points 
4>eff (—1, 0, 0) and (1, 0, 0) yields to a critical Jacobi con¬ 
stant Cl , which can be used to define a dimensionless energy 
parameter as 


^ Cl-E 
Cl 


(16) 


where E < 0 is some other value of the Jacobian. The di¬ 
mensionless energy parameter C makes the reference to en¬ 
ergy levels more convenient. For E — Cl the equipoten- 
tial surface encloses the critical volume, while for a Jaco¬ 
bian value E > Cl, or in other words C > 0, the equipo- 
tential surface is open and consequently stars can escape 
from the cluster. In the case of the isochrone potential 
Cli — —3.543290584540038, while for the Hernquist poten¬ 
tial we have Cl u — —3.537211521390788. 

It would be very informative to study and compare the 


mass density derived from the potentials we use to describe 
the star cluster. Since all potentials are spherically symmet¬ 
ric the mass density can be easily computed using the fol¬ 
lowing version of the Poisson’s equation 


p(r) = J-1A(V 

; 4t tG r 2 dr \ dr J 


(17) 


After trivial computations we have 


Pp(r) 



5/2 


pi m = 


3M c i ((n + a)a 2 - r 2 (n + 3a)) 
47T (n + a) 3 a 3 


Pa(r) = 


Mci r H 1 
27T r (r + r H ) 3 ’ 


(18) 


where a = (r 2 + r 2 ) 1/<2 , while pp(r), pi(r), PhM are the 
densities corresponding to the Plummer, isochrone and 
Hernquist potential, respectively. In Fig. 1 we present the 
evolution of the mass density of the three potentials as a 
function of the radius r. For the Plummer model we take 
M c i = 2.7 and cp = 0.38. The value of cp was calculated as 
the other two scale radii. It is seen that for distances very 
close to the center (r < 0.1) the density of the Plummer and 
isochrone potentials remain constant having finite-density 
core, while on the other hand the density of the Hernquist 
model increases due to the gentle power-law cusp in the 
center. For intermediate distances (0.1 < r < 1) the evolu¬ 
tion of all three density profiles is almost the same. We also 
observe that for large radii the asymptotic behavior of the 
isochrone and the Hernquist mass density profiles is identical 
varying like 1/r 4 , while the density of the Plummer poten¬ 
tial drops more quickly following the power law 1/r 5 . Before 
closing this section we would like to add that Miocchi et al. 
(2013) showed that King (King 1966) and Wilson (Wilson 
1975) density profiles fit globular clusters well. Both King 
and Wilson models have, in contrast to the Hernquist model, 
finite densities in the core. 


3 COMPUTATIONAL METHODS 

In order to investigate the orbital content of the star clus¬ 
ter, we need to define samples of initial conditions of orbits 
whose properties (bounded or escaping motion) will be iden¬ 
tified. For this purpose, we define for each value of the total 
orbital energy (all tested energy levels are above the escape 
energy), dense uniform grids of 1024 x 1024 initial condi¬ 
tions regularly distributed in the area allowed by the value 
of the energy. Our investigation takes place both in the con¬ 
figuration (x,y) and in the phase (x,x) space for a better 
understanding of the escape mechanism. Furthermore, the 
grids of initial conditions of orbits whose properties will be 
explored are defined as follows: For the configuration (x, y) 
space we consider orbits with initial conditions (xo,2/o) with 
yo = 0, while the initial value of x’o is always obtained from 
the Jacobian integral (15) as x'o = x(xo, yo , y‘o, E) > 0. Simi¬ 
larly, for the phase (x, x) space we consider orbits with initial 
conditions (xo,£o) with yo — 0, while this time the value of 
yo is obtained from the Jacobian integral (15). 
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The equations of motion as well as the variational equa¬ 
tions for the initial conditions of all orbits were integrated 
using a double precision Bulirsch-Stoer FORTRAN 77 algo¬ 
rithm (e.g., Press et ah 1992) with a small time step of or¬ 
der of 10 -2 , which is sufficient enough for the desired accu¬ 
racy of our computations (i.e., our results practically do not 
change by halving the time step). Here we should empha¬ 
size, that our previous numerical experience suggests that 
the Bulirsch-Stoer integrator is both faster and more ac¬ 
curate than a double precision Runge-Kutta-Fehlberg algo¬ 
rithm of order 7 with Cash-Karp coefficients (e.g., Darriba 
2012). Throughout all our computations, the Jacobian en¬ 
ergy integral (Eq. (15)) was conserved better than one part 
in 10 —11 , although for most orbits it was better than one 
part in 10 -12 . 

An issue of paramount importance in dynamical sys¬ 
tems with escapes is the determination of the position as 
well as the time at which an orbit escapes from the system. 
For all energy values smaller than the critical value Cl (es¬ 
cape energy), the 3D equipotential surface is closed. On the 
other hand, when E > Cl the equipotential surface is open 
and extend to infinity. The value of the energy itself how¬ 
ever, does not furnish a sufficient condition for escape. An 
open equipotential surface consists of two branches forming 
channels through which an orbit can escape to infinity. It 
was proved (Churchill et al. 1979) that in Hamiltonian sys¬ 
tems at every opening there is a highly unstable periodic 
orbit close to the line of maximum potential which is called 
a Lyapunov orbit. Such an orbit reaches the Zero Veloc¬ 
ity Curve (ZVC), on both sides of the opening and returns 
along the same path thus, connecting two opposite branches 
of the ZVC. Lyapunov orbits are very important for the es¬ 
capes from the system, since if an orbit intersects any one of 
these orbits with velocity pointing outwards moves always 
outwards and eventually escapes from the system without 
any further intersections with the surface of section (e.g., 
Contopoulos 1990). When E — Cl the Lagrangian points 
exist precisely but at E > Cl an unstable Lyapunov peri¬ 
odic orbit is located close to each of these two points (e.g., 
Henon 1969). The Lagrangian points are saddle points of 
the effective potential, so when E > Cl, a star must pass 
close enough to one of these points in order to escape. Thus, 
in the case of a tidally limited star cluster the escape crite¬ 
rion is purely geometric. In particular, escapers are defined 
to be those stars moving in orbits beyond the tidal radius 
thus passing one of the two Lagrangian points with velocity 
pointing outwards (Fukushige V Heggie 2000). 

The configuration and the phase space are divided into 
the bounded and unbounded regions. Usually, the vast ma¬ 
jority of the bounded space is occupied by initial conditions 
of regular orbits forming stability islands. In many systems 
however, trapped chaotic orbits have also been observed. 
Therefore, we decided to distinguish between regular and 
chaotic trapped orbits. Over the years, several chaos indica¬ 
tors have been developed in order to determine the character 
of orbits. In our case, we choose to use the Smaller ALing- 
ment Index (SALI) method. The SALI (Skokos 2001) has 
been proved a very fast, reliable and effective tool, which is 
defined as 

SALI(t) = min(d_,d+), (19) 

where d- = ||wi(£) — W 2 (f) || and = ||wi(£) + W 2 (t)|| are 


the alignments indices, while wi(f) and W 2 (f), are two devi¬ 
ation vectors which initially point in two random directions. 
For distinguishing between ordered and chaotic motion, all 
we have to do is to compute the SALI along time interval 
tmax of numerical integration. In particular, we track simul¬ 
taneously the time-evolution of the main orbit itself as well 
as the two deviation vectors wi(f) and W 2 (t) in order to 
compute the SALI. The time-evolution of SALI strongly de¬ 
pends on the nature of the computed orbit since when an 
orbit is regular the SALI exhibits small fluctuations around 
non zero values, while on the other hand, in the case of 
chaotic orbits the SALI after a small transient period it 
tends exponentially to zero approaching the limit of the ac¬ 
curacy of the computer (10 -16 ). Therefore, the particular 
time-evolution of the SALI allow us to distinguish fast and 
safely between regular and chaotic motion. Nevertheless, we 
have to define a specific numerical threshold value for deter¬ 
mining the transition from order to chaos. After conducting 
extensive numerical experiments, integrating many sets of 
orbits, we conclude that a safe threshold value for the SALI 
is the value 10 -8 . In order to decide whether an orbit is reg¬ 
ular or chaotic, one may follow the usual method according 
to which we check after a certain and predefined time inter¬ 
val of numerical integration, if the value of SALI has become 
less than the established threshold value. Therefore, if SALI 
^ 10 -8 the orbit is chaotic, while if SALI > 10 -8 the orbit 
is regular thus making the distinction between regular and 
chaotic motion clear and beyond any doubt. For the com¬ 
putation of SALI we used the LP-VI code (Carpintero et al. 
2014), a fully operational routine which efficiently computes 
a suite of many chaos indicators for dynamical systems in 
any number of dimensions. 

In our calculations, we set 10 4 time units as a maximum 
time of the numerical integration. The vast majority of or¬ 
bits (regular and chaotic) however, need considerable less 
time to find one of the two exits in the limiting surface and 
eventually escape from the system (obviously, the numerical 
integration is effectively ended when an orbit passes through 
one of the escape channels and escapes). Nevertheless, we de¬ 
cided to use such a vast integration time just to be sure that 
all orbits have enough time in order to escape. Remember, 
that there are the so called “sticky orbits” which behave as 
regular ones during long periods of time (see also Chapter 47 
in Contopoulos (2004)). The SALI method can easily iden¬ 
tify sticky orbits. In fact sticky orbits are those which hold 
intermediate values of SALI (10 -4 < SALI < 10 -8 ) for long 
time intervals (see e.g., Skokos et al. (2004)). It should be 
clarified, that orbits which do not escape after a numerical 
integration of 10 4 time units are considered as non-escaping 
or trapped regarding the value of SALI. In fact, orbits with 
escape periods larger than 10 4 time units are completely ir¬ 
relevant to our investigation since they lack physical mean¬ 
ing. 


4 NUMERICAL RESULTS 

In our investigation we seek to determine which orbits es¬ 
cape and which remain trapped, distinguishing simultane- 
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ously between regular and chaotic bounded motion 1 . More¬ 
over, two additional aspects of the orbits will be examined: 
(i) the channel or exit through which a star escapes and (ii) 
the time-scale of the escape (we shall also use the terms es¬ 
cape period or escape rate). In particular, we shall explore 
these dynamical quantities for various values of the total 
orbital energy, always within the interval C G [0.001,0.1]. 
As in Paper I, we shall deal only with the two-dimensional 
(2D) case of the star cluster model, therefore everywhere 
z = z = 0. 

Our preliminary analysis suggests that apart from the 
escaping orbits there is also a considerable amount of 
bounded orbits. In general terms, the majority of non¬ 
escaping bounded regions corresponds to initial conditions 
of regular orbits, where the adelphic integral of motion is 
present, restricting their accessible phase space and there¬ 
fore hinders their escape. However, there are also chaotic 
orbits which do not escape within the predefined interval 
of 10 4 time units and remain trapped for vast periods un¬ 
til they eventually escape to infinity. Therefore, we decided 
to classify the initial conditions of orbits in both the con¬ 
figuration and the phase space into four main categories: 
(i) orbits that escape through Li, (ii) orbits that escape 
through L 2 , (iii) non-escaping bounded regular orbits and 
(iv) trapped chaotic orbits. Additional numerical computa¬ 
tions reveal that the non-escaping regular orbits are mainly 
loop orbits 2 * for which the adelphic integral applies, while 
other types of secondary resonant orbits are also present. 
In Fig. 2a we present an example of a 3:3 resonant orbit, 
while in Fig. 2b a characteristic example of a secondary 1:3 
resonant orbit is shown. The n : m notation we use for the 
regular orbits is according to Carpintero & Aguilar (1998) 
and Zotos & Carpintero (2013), where the ratio of those in¬ 
tegers corresponds to the ratio of the main frequencies of 
the orbit, where main frequency is the frequency of great¬ 
est amplitude in each coordinate. Main amplitudes, when 
having a rational ratio, define the resonances of an orbit. 
Finally in Figs. 2c-d we observe two orbits escaping through 
L\ (exit channel 1) and L 2 (exit channel 2), respectively. 
The unstable Lyapunov orbits near the Lagrangian points 
are shown in red. Both regular orbits shown in Figs. 2a-b are 
computed until t = 100 time units, while on the other hand 
the escaping orbits presented in Figs. 2c-d were calculated 
for 5 time units more than the corresponding escape period 
in order to visualize better the escape trail. All orbits have 
yo = xo = 0, while the value of yo is obtained from the Ja¬ 
cobian integral where C = 0.01. In Table 1 we provide the 
type, the exact initial condition xo and the escape period 
for all the depicted orbits. For modelling the star cluster 
the isochrone potential was used, while for the case of the 
Hernquist potential the types of orbits are similar. In Fig. 2 
we observe the two openings (exit channels) through which 

1 Generally, any dynamical method requires a sufficient time in¬ 
terval of numerical integration in order to distinguish safely be¬ 
tween ordered and chaotic motion. Therefore, if the escape period 
of an orbit is very low or even worse if the orbit escapes directly 
from the system then, any chaos indicator will fail to work prop¬ 
erly due to insufficient integration time. Hence, it is pointless to 
speak of regular or chaotic escaping orbits. 

2 Orbits which form a rosette or a simple circular closed path are 

known as loop orbits (see also Chapter 2 in Contopoulos (2002)). 


Table 1 . Type, initial conditions, and integration time (tint) of 
the orbits shown in Fig. 2(a-d). 


Figure 

Type 

xo 

tint 

2a 

3:3 

-0.625 

100 

2b 

1:3 

-0.911 

100 

2c 

chaotic 

0.620 

76 

2d 

chaotic 

0.670 

79 


the particles can leak out. In fact, we may say that these 
two exits act as hoses connecting the interior region of the 
star cluster where —rt<x<rt with the “outside world” 
of the exterior region. Exit channel 1 (negative ^-direction) 
indicates escape towards the galactic center, while channel 2 
(positive ^-direction) indicates escape towards infinity. The 
forbidden regions of motion within the banana-shaped iso¬ 
lines are shown in the same figure with gray. 

A simple qualitative way for distinguishing between 
regular and chaotic motion in Hamiltonian systems is by 
plotting the successive intersections of the orbits using a 
Poincare Surface of Section (PSS) (Henon and Heiles 1964). 
In particular, a PSS is a two-dimensional (2D) slice of the 
entire four-dimensional (4D) phase space. The following Fig. 
3(a-d) shows the PSSs for both potentials at the critical Ja¬ 
cobi energy Cl- Fig. 3a corresponds to the configuration 
(x,y) space showing orbits crossing y — 0 with x > 0, while 
orbits crossing y — 0 with y > 0 at the phase (x,x) space 
are presented in Fig. 3b for the isochrone potential. Things 
are similar in Figs. 3(c-d) for the Hernquist potential. In 
both cases, the initial conditions of the orbits are integrated 
forwards in time plotting a dot at each crossing through the 
surface of section. The outermost black thick curve circum¬ 
scribing each PSS is the limiting curve which in the configu¬ 
ration space is defined as (x, y, z = 0) = E, while in the 
phase space is given by 

f(x,x) = ii 2 + $eff (x,y = 0,z = 0) = E. (20) 

It is seen, that the vast majority of the surfaces of section is 
covered by a unified and extended chaotic sea however, sta¬ 
bility regions with initial conditions corresponding to regular 
orbits are also present. In these stability islands the quasi- 
periodic orbits admit a third integral of motion, also known 
as “adelphic integral” (see e.g., Contopoulos (1963)), which 
hinders the test particles (stars) from escaping. Moreover, 
we observe that some areas in the chaotic sea of the (x, y) 
plane are less densely occupied than others. We remark that 
in the following grid exploration in the configuration and 
phase space, we do not consider all initial conditions in the 
same chaotic domain as one the same orbit. 

4.1 The configuration (x,y) space 

Our exploration begins in the configuration (x, y) space. Fig. 
4 shows the orbital structure of the configuration space for 
three energy levels, where the four basic types of orbits are 
indicated with different colors. Specifically, gray color cor¬ 
responds to regular non-escaping orbits, white color corre¬ 
sponds to trapped chaotic orbits, green color corresponds to 
orbits escaping through channel 1, while the initial condi¬ 
tions of orbits escaping through exit channel 2 are marked 
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Figure 2. Characteristic examples of orbits in the star cluster where the isochrone potential is used, (a): 3:3 resonant orbit, (b): 1:3 
resonant orbit, (c): orbit escaping through L i, (d): orbit escaping through L 2 . The red lines near the Lagrangian points denote the 
unstable Lyapunov orbits. 


with red color. In the left column we present the results 
where the isochrone potential is used, while in the right col¬ 
umn we see the outcomes of the Hernquist potential. Our 
results will be compared with the corresponding ones of the 
Plummer model shown in the right column of Fig. 5 of Pa¬ 
per I 3 . We observe that for all energy levels there are many 
similarities between the three models but there are also sig¬ 
nificant differences. In particular, it is seen that the main 
stability islands containing the initial conditions of all the 
1:1 resonant orbits are almost the same everywhere. How¬ 
ever the structure of the secondary stability islands is very 

3 In Paper I the total mass of the star cluster was M c \ = 2.2, 
while in the cases of the isochrone and the Hernquist potential 
the corresponding value is M c \ = 2.7. However, our calculations 
suggest that this minor deviation in the total mass of the star 
cluster does not significantly affect the orbital structure. There¬ 
fore, we decided not to reconstruct the orbital diagrams of the 
Plummer case but instead use those of Paper I. 


different. Specifically one can see that in both the isochrone 
and the Hernquist potentials the set of four small islands 
of the Plummer case corresponding to the 4:4 resonant or¬ 
bits which are located around the lower 1:1 region is now 
absent. Moreover, the set of the 3:3 stability islands around 
the upper 1:1 region is absent in the case of the Hernquist po¬ 
tential. Undoubtedly, the most striking difference concerns 
the 1:3 stability regions. In particular, in the isochrone case 
their domains are more extended than in the case of the 
Plummer potential. The 1:3 stability regions are much more 
extended in the Hernquist case where we observe that the 
three islands are almost connected to each other. With a 
more closer look we can distinguish a thin layer of initial 
conditions of trapped chaotic orbits separating the 1:3 re¬ 
gions. When C — 0.001 in all thee cases the vast majority 
of the ( x , y ) plane is highly fractal, while for higher energy 
levels the fractality of the plane is reduced and well-formed 
basins of escape dominate. When C — 0.1 (the highest en¬ 
ergy level studied) the 1:3 stability islands have disappeared 
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Figure 3. Poincare Surfaces of Section (PSS) at the critical Jacobi energy Cl for (first row): the isochrone potential and (second row): 
the Hernquist potential, (a and c left): The configuration (x,y) space where orbits crossing y = 0 with x > 0 and (b and d right): the 
phase ( x , x) space where orbits crossing y = 0 with y > 0. All types of regular orbits are identified. 


and the vast majority of the configuration plane is occupied 
by several basins of escape. In all three cases the structure of 
the escape basins is quite similar and only small differences 
can be mentioned (i.e, the fractal structures at the lower 
part (y < 0) of the configuration plane are more prominent 
in the Hernquist case when C = 0.1). 

The first column of Fig. 5 shows the the time-evolution 
of the x and y coordinates of a fast escaping orbit in the 


Hernquist potential with initial conditions: xq — —0.52, 
yo = 0, yo = 0, while the value of xo is obtained from 
the Jacobian integral (15) when C — 0.001. The horizon¬ 
tal dashed red lines in the upper panel denote the position 
of the two Lagrangian points delimiting the transition from 
trapped to escaping motion. In the same figure we see the 
time-evolution of SALI as well as the time-evolution of the 
error in the Jacobian integral (A E — (\E — Eq\)/Eq). We 
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Figure 4. Orbital structure of the configuration ( x , y) space. Left column: case of the isochrone potential; Right column: case of the 
Hernquist potential. Top row: C — 0.001; Middle row: C — 0.01; Bottom row: C — 0.1. The green regions correspond to initial conditions 
of orbits where the stars escape through L i, red regions denote initial conditions of orbits where the stars escape through L 2 , gray areas 
represent stability islands of regular bounded orbits, while initial conditions of trapped chaotic orbits are shown in white. 
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Figure 5. Time-evolution of x and y coordinates, the SALI and the integration error in the Jacobian integral of (first column): a fast 
escaping orbit; (second column): a sticky orbit; (third column): a trapped chaotic orbit in the case of the Hernquist potential. 


see that the orbit need only about 110 time units in order to 
cross the SALI threshold value, while it escapes from channel 
2 after about 246 time units. The error in the Jacobian inte¬ 
gral proves the accuracy of our numerical integrator. In the 
second column of Fig. 5 we present the time-evolution of a 
stick orbit with initial conditions: xo = 0, yo = 0.062, y 0 = 0, 
while the value of xo is obtained from the Jacobian integral 
(15) when C — 0.001. According to the time-evolution of 
SALI this orbit behave as regular (SALI > 10 _4 ) for the first 
about 600 time units. Then in the interval 600 < t < 900 it 
displays intermediate values of SALI (sticky period) before 
it escapes form channel 1 after about 2020 time units. 

In the case of the Hernquist potential (see right column 
of Fig. 4) we identified some trapped chaotic orbits. At this 
point we would like to try to interpret and also justify the 
phenomenon of trapped chaos. By the term “trapped chaos” 
we refer to the case where a chaotic orbit remains trapped 
for a vast time interval inside an open equipotential surface. 
It should be emphasized and clarified that these trapped 
chaotic orbits cannot be considered by no means neither as 
sticky orbits nor as super sticky orbits 4 with sticky periods 
larger than 10 4 time units. Sticky orbits are those who be¬ 
have regularly for long time periods before their true chaotic 


4 A super sticky orbit is an orbit with a sticky period larger than 
10 4 time units. The sticky period is defined as the time interval 
in which SALI has intermediate values (10 -4 < SALI < 10 -8 ) 
before it crosses the threshold value of 10 -8 . 


nature is fully revealed. In our case on the other hand, this 
type of orbits exhibit chaoticity very quickly as it takes no 
more than about 100 time units for the SALI to cross the 
threshold value (SALI <C 10 -8 ), thus identifying beyond any 
doubt their chaotic character. Chaotic orbits do not admit a 
third integral of motion, therefore according to the classical 
chaos theory for time t —> oo they will fill all the available 
space inside the equipotential surface, unless they escape. 
In other words, given enough time of numerical integration 
all trapped chaotic orbits will eventually pass through an 
exit and escape from the star cluster. To prove this, we 
chose such an initial condition and let the time running. 
In the third column of Fig. 5 we present the time-evolution 
of the x and y coordinates of an orbit with initial condi¬ 
tions: xo = —0.37582418, yo = —0.10329670, yo — 0, while 
the value of xo is obtained from the Jacobian integral (15) 
when C — 0.01. This particular orbit is trapped chaotic and 
as we can see a huge time of numerical integration, about 
558000 time units is required so that the orbit can escape 
from channel 2. It should be emphasized that for this or¬ 
bit the SALI crosses the threshold value after about 150 
time units and its values remain in the chaotic zone (SALI 
i 10 -8 ) during the entire time interval before the escape. 
We observe in the last panel that the error in the Jacobian 
integral grows during the vast numerical integration but al¬ 
ways in reasonable values. We experimented with several 
trapped chaotic orbits and the obtained results were very 
similar. Stars moving in chaotic orbits restricted inside the 
interior region of the cluster having large escape times have 
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Figure 6. Distribution of the escape times t esc of the orbits on the configuration (x, y) space. Left column: case of the isochrone potential; 
Right column: case of the Hernquist potential. Top row: C = 0.001; Middle row: C = 0.01; Bottom row: C = 0.1. The darker the color, 
the larger the escape time. Initial conditions of non-escaping regular orbits and trapped chaotic orbits are shown in white. 


also been reported in previous investigations (see e.g., Fig. 
11 in Fukushige & Heggie 2000). 

Perhaps the presence of a small percentage of trapped 
chaotic orbits in the case of the Hernquist potential is related 
to the presence of the cusp near the center of the potential 
where the density of Eq. (18) gets singular and the velocities 
become arbitrarily high. Additional numerical calculations 
indicate that the phenomenon of trapped chaos is not an 
artifact of the numerical integration since the areas covered 
by trapped chaotic orbits remain unperturbed when we use 
smaller time step (0.005 or 0.001) in the numerical integra¬ 
tion. One may suspect that the trapped chaotic orbits might 
pass close to the origin at (x,y) = (0, 0), where the singular 
cusp is located and they might escape relatively fast because 
of an energy error in the Jacobian since the velocities can 


become relatively large near the singularity. Our numerical 
calculations however, indicate that this scenario is not plau¬ 
sible. These orbits do eventually escape but only after vast 
time intervals of more than 10 5 time units. 

Another possible reason for the phenomenon of trapped 
chaos observed in the Hernquist model might be the force 
acting on stars near the center of the star cluster. For both 
the Plummer and the isochrone potentials the force near the 
center vanishes, while in the case of the Hernquist potential 
on the other hand it is non zero. In particular, according 
Hernquist (1990) 


F(r) = 


dr 


GMci 

('r + m ) 2 


( 21 ) 


For r — 0 we have that F{r — 0) — —232.242. 
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Figure 7. Evolution of the percentages of escaping and non-escaping orbits when varying the energy parameter C on the configuration 
(x,y) space for the case of the (a-left): isochrone potential and (b-right): Hernquist potential. 


In the following Fig. 6 we present how the escape times 
t esc of orbits are distributed on the configuration ( x , y ) space 
for both potentials. Light reddish colors correspond to fast 
escaping orbits with short escape periods, dark blue/purple 
colors indicate large escape rates, while white color denote 
both non-escaping regular and trapped chaotic orbits. It 
is evident, that orbits with initial conditions close to the 
boundaries between the escape basins, that is, the fractal 
areas of the plots need significant amount of time in order 
to escape from the cluster, while on the other hand, inside 
the basins of escape where there is no dependence on the 
initial conditions whatsoever, we measured the shortest es¬ 
cape rates of the orbits. We observe that for C = 0.001 the 
escape periods of orbits with initial conditions in the fractal 
basin boundaries are huge corresponding to tens of thou¬ 
sands of time units. This phenomenon is anticipated because 
in this case the width of the escape channels is very small 
and therefore, the orbits should spend much time inside the 
equipotential surface until they find one of the two openings 
and eventually escape to infinity. One may observe, that as 
the value of the energy increases however, the escape chan¬ 
nels become more and more wide, while the escape times of 
all orbits (inside the exit basins and in the fractal regions) is 
reduced considerably. Our calculations suggest that statisti¬ 
cally the escape times in the case of the Hernquist potential 
are lower (almost half) than in the case where the star clus¬ 
ter is modelled by the isochrone potential. Moreover, we see 
that for C — 0.01 and C — 0.1 the escaping orbits with ini¬ 
tial conditions near the boundaries of the stability islands 
have high escape times (more than 1000 time units). This is 
true because near the boundaries of the stability islands is 
usually the location of sticky orbits which need a sufficient 
time interval before they reveal their true chaotic nature and 
finally escape from the system (see second column of Fig. 5). 

The evolution of the percentages of all types of orbits 
on the configuration space for both potentials when the di¬ 
mensionless energy parameter C varies is presented in Fig. 
7(a-b). Due to an unreasonably difficulty in the data analysis 


we decided to investigate in the configuration space (similar 
approach in the following phase space) only five different val¬ 
ues of C . The case of the isochrone potential is shown in Fig. 
7a. When C — 0.001, that is an energy level just above the 
escape energy, about one fourth of the total plane is covered 
by initial conditions of non-escaping regular orbits, while the 
rest 75% is devoted to escaping orbits. In particular the rates 
of both channels are identical suggesting that in this case the 
configuration plane is highly fractal. As we proceed to higher 
energy levels the percentage of the non-escaping regular or¬ 
bits is reduced and when C — 0.1 is only about 12%. On 
the contrary we observe that always the escaping orbits is 
the most populated family. With increasing energy the rates 
of both channels (exits) start to diverge and especially for 
C > 0.01 the percentage of escapers through L\ decreases, 
while that of escaping orbits through L 2 increases linearly. 
At the highest energy level studied (C = 0.1), escaping or¬ 
bits through channel 1 correspond to about 18% of the ( x , y ) 
plane, while escapers through channel 2 dominate covering 
the rest 70%. In the case of the Hernquist potential (Fig. 7b) 
the profile of the curves is similar however there are some 
important differences: (i) in the interval 0.001 < C < 0.01 
non-escaping regular orbits is the most populated type of 
orbits; (ii) the percentage of escaping orbits trough channel 
1 exhibits only a minor decrease from 28% to 24%; (iii) es¬ 
caping orbit through L 2 dominate only for C > 0.01. At this 
point we should point out that in both cases trapped chaotic 
orbits possess always a very weak percentage (less than 1%) 
and this is the main reason why this evolution of this type is 
not included in the diagrams. Thus taking into account all 
the above-mentioned analysis we may reasonably conclude 
that in the Hernquist models there is a greater percentage 
of bounded regular orbits mainly due to the strong presence 
of the 1:3 family. Furthermore, in both cases for relatively 
high energy levels escaping orbits through exit 2 dominate. 
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Figure 8. Orbital structure of the phase (x, x) space. Left column: case of the isochrone potential; Right column: case of the Hernquist 
potential. Top row: C — 0.001; Middle row: C = 0.01; Bottom row: C — 0.1. The green regions correspond to initial conditions of orbits 
where the stars escape through L i, red regions denote initial conditions of orbits where the stars escape through L 2 , gray areas represent 
stability islands of regular bounded orbits, while initial conditions of trapped chaotic orbits are shown in white. 
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Figure 9. Distribution of the escape times t esc of the orbits on the phase (x,x) space. Left column: case of the isochrone potential; 
Right column: case of the Hernquist potential. Top row: C = 0.001; Middle row: C = 0.01; Bottom row: C = 0.1. The darker the color, 
the larger the escape time. Initial conditions of non-escaping regular orbits and trapped chaotic orbits are shown in white. 


4.2 The phase (x,x) space 

Our investigation continues in the phase (x,x) space where 
we follow the same numerical approach as discussed previ¬ 
ously. In Fig. 8 we depict the orbital structure of the (x,x) 
plane for the isochrone and Hernquist cases and for the three 
energy levels, using different colors in order to distinguish 
between the four main types of orbits (non-escaping regular; 
trapped chaotic; escaping through L\ and escaping through 
L 2 ). This time our results will be compared with the cor¬ 
responding ones of the Plummer model shown in the left 
column of Fig. 5 of Paper I. The observed regular areas cor¬ 
respond mainly to retrograde 1:1 resonant orbits (see also 
Fukushige & Heggie 2000; Ernst et al. 2007) (i.e., when a 
star revolves clockwise around the cluster in the opposite 


sense with respect to the motion of the cluster around the 
parent galaxy), while there is also a smaller stability island 
of prograde counterclockwise 1:1 resonant orbits in the right 
part {x > 0) of the plane. As for the other resonant orbits 
the similarities and differences are the same with that ex¬ 
plained in the previous subsection. However there is a new 
significant difference. In particular, we observe that in the 
case of the Hernquist potential stars can obtain much greater 
velocity near the center due to the power-law cusp of the po¬ 
tential. Regarding the structure of the fractal areas and the 
several escape basins we may say that in general terms there 
are no important differences between the three cases. 

Here we must point out, that the phase (x, x) plane 
shown in Fig. 8 is not a classical PSS, simply because escap¬ 
ing orbits in general, do not intersect the y — 0 axis after 
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Figure 10. Evolution of the percentages of escaping and non-escaping orbits when varying the energy parameter C on the phase ( x , x) 
space for the case of the (a-left): isochrone potential and (b-right): Hernquist potential. 


a certain time, thus preventing us from defying a recurrent 
time. A classical PSS exists only if orbits intersect an axis, 
like y — 0, at least once within a certain time interval. Nev¬ 
ertheless, in the case of escaping orbits we can still define 
local surfaces of section which help us to understand the 
orbital behavior of the dynamical system. 

We see in the ( x , x) surface of section shown in Fig. 3d 
for the Hernquist model that the invariant curves around 
the 1:3 periodic orbit located on the x axis block the exit 
close to L\. Of course in this case the channel is closed due 
to the fact that the surface of section is at the critical value 
of the Jacobi integral. However when E > Cl the stars mov¬ 
ing in trapped chaotic orbit, after more than 10 5 time units, 
escape between these invariant curves and the boundary of 
the surface of section on the top or bottom side and, after 
than, they reach the exist at L\. Our computations suggest 
that trapped chaotic orbits exist only for C = 0.001 and 
C — 0.01 when the stability islands of the 1:3 resonance 
are present. On the other hand, for C — 0.1 where the 1:3 
stability islands disappear we did not identify any initial 
conditions of trapped chaotic orbits. Therefore we may say 
that the presence of the 1:3 stability island may be respon¬ 
sible, in a way, for the observed trapped chaos. It should be 
noted however, that the 1:3 stability island is also present in 
the case of the isochrone model, where no trapped chaotic 
orbits were reported. 

The distribution of the escape times t esc of orbits on 
the phase space for both cases is shown in Fig. 9. One may 
observe that the results are very similar to those presented 
earlier in Fig. 6, where we found that orbits with initial 
conditions inside the escape basins have the smallest escape 
rates, while on the other hand, the longest escape times 
correspond to orbits with initial conditions in the fractal 
regions of the plots or in the boundaries of the stability 
islands. At this point, we would like to point out that the 
basins of escape can be easily distinguished being the regions 
with intermediate greenish colors indicating fast escaping 
orbits. Indeed, our numerical calculations suggest that orbits 


with initial conditions inside the basins need no more that 
10 time units in order to escape from the star cluster. 

Fig. 10(a-b) shows how the percentages of all types of 
orbits on the phase space evolve for both cases when the 
energy parameter varies in the interval C G [0.001, 0.1]. For 
the isochrone potential we observe in Fig. 10a that for low 
energy levels non-escaping regular orbits occupy about 36% 
of the phase space. As the value of the energy increases 
however, the portion of the ordered orbits is reduced and 
when C — 0.1 it corresponds to only 24% of the (x, x) plane. 
It is also seen that the rates of the escaping orbits evolve 
almost identically and only at relatively high values of the 
energy the rates slightly diverge, being the escapers through 
exit 2 more populated. Things are quite similar in the case 
of the Hernquist potential. Indeed in Fig. 10b we see that in 
the interval 0.001 < C < 0.01 bounded regular orbits cover 
about 46.5% of the phase space, while for larger values of 
the energy there is a steep fall of the rate to about 23% when 
C — 0.1. Furthermore, the rates of the escaping orbits are 
exactly the same everywhere. Initially escaping orbits share 
about 52% of the phase space, while for C > 0.01 their rates 
exhibit a significant growth and at the highest energy level 
studied we found that 77% of the phase space is covered 
by initial conditions of escaping orbits. Once again, as we 
found in the configuration space, the rate of trapped chaotic 
orbits is extremely low (always less than 1%). Therefore one 
can argue that in the phase space both exit channels are 
equiprobable throughout the examined energy interval. 

4.3 An overview analysis 

The color-coded grids that discussed earlier in the config¬ 
uration (x,y) as well as in the phase (x,x) space provide 
sufficient information on the phase space mixing for only 
a fixed value of the Jacobi constant. Henon however, back 
in the late 60s (e.g., Henon 1969), introduced a new type 
of plane which can provide information not only about sta¬ 
bility and chaotic regions but also about areas of bounded 
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Figure 11. (first row): Orbital structure of the (x, C)-plane; (second row): the distribution of the corresponding escape times of the 2D 
orbits, (first column): Plummer potential; (second column): Isochrone potential; (third column): Hernquist potential. The color codes 
are the exactly same as in Figs. 4 and 6. 
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Figure 12. Evolution of the percentages of escaping and non-escaping orbits on the (x,C )-plane as a function of the dimensionless 
energy parameter C for the (a): Plummer potential; (b): isochrone potential; (c): Hernquist potential. 


and escaping orbits using the section y = x = 0, y > 0. In 
other words, all the 2D orbits of the stars of the cluster are 
launched from the x-axis with x = xo, parallel to the y -axis 
(■y — 0). Consequently, in contrast to the previously pre¬ 
sented types of planes, only orbits with pericenters on the 
x-axis are included and therefore, the value of the dimension¬ 
less energy parameter C can be used as an ordinate. In this 
way, we can monitor how the energy influences the overall 
orbital structure of the system using a continuous spectrum 


of energy values rather than few discrete energy levels. In 
the first row of Fig. 11 we present the orbital structure of 
the (x, C )-plane when C G [0.001,0.1], while in the second 
column of the same figure the distribution of the correspond¬ 
ing escape times of orbits is depicted. Unfortunately, Paper 
I does not contain such an analysis, so we had to compute 
also the case of the Plummer potential for M c \ — 2.7 and 
cp — 0.38. 

It is seen that in all cases for low values of the energy, 
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Figure 13. (a-left): The average escape time of orbits < t esc > and (b-right): the percentage of the (x,C) planes covered by the escape 
basins as a function of the dimensionless energy parameter C for all three cases. 


above but close enough to the critical escape energy, there 
is a highly fractal structure between the stability islands. 
As the value of the energy increases however, the fractality 
of the planes decreases and several basins of escape start 
to emerge, mainly at the outer parts of the planes. We can 
also identify three main stability islands. The one on the 
left side of the planes correspond to the 1:3 resonant orbits, 
the middle stability islands contains the 1:1 retrograde or¬ 
bits, while that on the right side of the plane corresponds 
to 1:1 prograde orbits. The most weak 1:3 stability region is 
observed in the Plummer case, while in the case of the Hern- 
quist potential we have the most prominent presence. In all 
three cases the 1:3 resonance seems to vanish for C > 0.03. 
The stability island of the 1:1 retrograde orbits is more or 
less the same in all three cases. On the other hand, in the 
case of the Hernquist potential the stability island of the 
1:1 prograde orbits seems to disappear for C > 0.1, while 
at the other two cases is still present. For the Plummer and 
the isochrone cases we can distinguish the two 3:3 stabil¬ 
ity regions that embrace the main 1:1 prograde island. It 
is evident from the results shown in Figs, ll(d-f) that the 
escape times of the orbits are strongly correlated to the es¬ 
cape basins. Moreover, one may conclude that the smallest 
escape periods correspond to orbits with initial conditions 
inside the escape basins, while orbits initiated in the frac¬ 
tal regions of the planes or near the boundaries of stability 
islands have the highest escape rates. In all three cases the 
escape times of orbits are significantly reduced with increas¬ 
ing energy. Combining all the numerical outcomes presented 
in Figs, ll(d-f) we may say that the key factor that deter¬ 
mines and controls the escape times of the orbits is the value 
of the orbital energy (the higher the energy level the shorter 
the escape rates), while the fractality of the basin bound¬ 


aries varies strongly both as a function of the energy and of 
the spatial variable. 

It would be very illuminating to monitor the evolution 
of the percentages of the different types of orbits on the 
(, x , C ) planes as a function of the dimensionless energy pa¬ 
rameter. For this purpose we introduce Fig. 12(a-c) which 
contains the results for all three cases. We see that for the 
Plummer and the isochrone cases the percentage of bounded 
regular orbits starts at about 55% and gradually reduces 
with increasing energy, while in the case of the Hernquist 
potential the rate of non-escaping regular orbits is about 
58%. In the Plummer case bounded motion is always the 
most populated type of motion, while in the other two cases 
escaping orbits prevail at high enough values of the energy. 
In the interval 0.001 < C < 0.02 the percentages of escap¬ 
ing orbits seem to fluctuate around 22% in all three cases, 
while for larger energy levels the rates start to diverge and 
orbits escaping through exit channel 2 dominate. It should 
be emphasized that the greater divergence is observed in the 
Hernquist case, where escapers form L 2 occupy about 45% 
of the (x, C) plane when C = 0.1. 

The evolution of the average value of the escape time 
< t esc > of orbits as a function of the dimensionless energy 
parameter is given in Fig. 13(a) for the (x, C ) plane. It is 
seen, that the profiles of the Plummer and the isochrone po¬ 
tential almost coincide, while that of the Hernquist potential 
has always lower values. In particular, for the Plummer and 
the isochrone cases and for low values of energy, just above 
the escape value, the average escape time of orbits is about 
800 time units (only 350 for the Hernquist case), however it 
reduces rapidly tending asymptotically to zero which refers 
to orbits that escape almost immediately from the system. 
We feel it is important to justify this behaviour of the escape 
time. As the value of the Jacobi constant increases the two 
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escape channels (which are of course symmetrical) become 
more and more wide and therefore, orbits need less and less 
time in order to find one of the openings in the open equipo- 
tential surface and eventually escape from the system. This 
geometrical feature explains why for low values of energy or¬ 
bits consume large time periods wandering inside the open 
equipotential surface until they eventually locate one of the 
two exits and escape to infinity. Finally Fig. 13(b) shows 
the evolution of the percentage of the total area ( A ) on the 
(x, C) plane corresponding to basins of escape, as a func¬ 
tion of the dimensionless energy parameter. It is seen that 
for low values of the energy in all three cases the planes are 
highly fractal (A is almost zero). However, as we proceed 
to higher energy levels the degree of fractalization reduces 
and the area corresponding to basins of escape start to grow 
rapidly. In the interval 0.001 < C < 0.03 A% is common 
in all three cases, while only for C > 0.03 the profile of the 
Hernquist potential start to diverge. Eventually, at very high 
energy levels (C — 0.1) the fractal domains are significantly 
confined and therefore the well formed basins of escape oc¬ 
cupy about 30% of the entire plane in the Plummer and the 
isochrone case and about 44% in the Hernquist case. Thus 
we may say that in the Hernquist model orbits escape more 
quickly and basins of escape are more extended that in the 
other two models. We assume that the deviations shown in 
Figs. 13(a-b) of the Hernquist potential with respect to the 
other two (Plummer and isochrone potentials) are mainly 
due to two reasons which are the only obvious differences 
between these three potentials: (i) the presence of the cen¬ 
tral cusp and (ii) the finite force near the center of the star 
cluster. 


5 DISCUSSION AND CONCLUSIONS 

In this work we tried to compare the orbital dynamics in star 
cluster models using three similar potentials for describing 
the properties of the spherically symmetric star cluster. The 
work was initiated in Paper I where a Plummer potential was 
applied, while in the present work we investigated the escape 
dynamics in the cases of the isochrone and the Hernquist 
potential. As in Paper I, we restricted our exploration in 
the two-dimensional model thus using for all orbits z — z — 
0. We managed to distinguish between ordered/chaotic and 
trapped/escaping orbits and we also located the basins of 
escape leading to different exit channels, finding correlations 
with the corresponding escape times of the 2D orbits. 

We defined for several values of the Jacobi integral, 
dense uniform grids of 1024 x 1024 initial conditions (#o, yo) 
and (xo,x'o) regularly distributed in the area allowed by the 
energy level on the configuration and phase space, respec¬ 
tively and then we identified regions of order/chaos and 
bound/escape. For the numerical integration of the orbits 
in each type of grid, we needed about between 7 hours and 
11 days of CPU time on a Pentium Dual-Core 2.2 GHz PC, 
depending on the escape rates of orbits in each case. For 
each initial condition, the maximum time of the numerical 
integration was set to be equal to 10 4 time units however, 
when a particle escaped the numerical integration was ef¬ 
fectively ended and proceeded to the next available initial 
condition. 

The present paper provides quantitative information re¬ 


garding the escape dynamics in star cluster models. The 
main numerical results of our research can be summarized 
as follows: 

(i) In the Plummer case two types of secondary reso¬ 
nances (the 3:3 and the 4:4 resonant families) surround the 
main 1:1 stability islands. In the isochrone case only the 3:3 
resonant family is present, while in the Hernquist models 
both secondary resonant families are absent. 

(ii) The 1:3 resonant family is common in all three cases 
however, in the Hernquist potential the areas on the con¬ 
figuration and on the phase space covered by 1:3 resonant 
orbits are much more extended with respect to the other two 
cases. 

(iii) At relatively high values of the energy the stability 
islands corresponding to the 1:3 resonant orbits disappear 
in all three cases. In addition, the stability island of the 
1:1 prograde orbits is confined (isochrone case) or disappear 
(Hernquist case). 

(iv) It was observed that the structure of the several es¬ 
cape basins (thin elongated bands or well formed broad re¬ 
gions) is very similar in all three cases, while all differences 
regarding the particular shapes are very minor. 

(v) Our calculations revealed, that the escape times of or¬ 
bits are directly linked to the basins of escape. In particular, 
inside the basins of escape as well as relatively away from 
the fractal domains, the shortest escape rates of the orbits 
had been measured. On the other hand, the longest escape 
periods correspond to initial conditions of orbits either near 
the boundaries between the escape basins or in the vicinity 
of the stability islands. 

(vi) In the Hernquist case we identified a non zero amount 
of trapped chaotic orbits which are mainly located inside 
the 1:3 stability islands. Additional numerical calculations 
revealed that these orbits do eventually escape to infinity 
but only after vast time intervals with no physical meaning. 

(vii) In the case of the Hernquist models we measured 
the smallest average escape times in all tested energy lev¬ 
els, while in the cases of the Plummer and the isochrone 
potentials the average escape rates of orbits are larger. Fur¬ 
thermore, the area of the basins of escape is larger in the 
Hernquist models. We assume that this divergence is due to 
the presence of the central cusp in the Hernquist potential. 

Judging by the detailed and novel outcomes we may 
say that our task has been successfully completed. We hope 
that the present comparison analysis and the correspond¬ 
ing numerical results to be useful in the active field of the 
dissolution of tidally limited star clusters. 
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